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Abstract 

For Hill's equation on [0, co) we prove new characterizations of the spectral function p(A) and the 
spectral density function /(A) based on analysis involving a companion system of first order differ- 
ential equations as in [SJ [7]. A numerical algorithm is derived and implemented based on coefficient 
approximation. Results for several examples, including the Mathieu equation, are presented. 



1 Introduction 

In this paper we consider Hill's equation, henceforth referred to as the SL- equation, 

— y" + q{x)y = \y, < x < oo, 
where q{x) is real valued and periodic with period I, and we impose a boundary condition 

y(0) cosa + y'(0) sina = 

for some a € [0, n). 

Let a fundamental system of solutions for be defined for all A G (C by 



0(0, A) #),A) 
6>'(0,A) 0'(O,A) 



cos a — sin a 
sin a cos a 



we can define a unique Titchmarsh-Weyl to— function by (Im A =/= 0) 

9(x, A) + m(\)(j)(x, A) € L 2 (0,oo). 



(1.1) 
(1.2) 

(1.3) 
(1.4) 
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The spectral function is then defined for A e [A, oo) by the Titchmarsh-Kodaira formula 

1 f x 

pfA) = lira — / Im m(/z + ie) du, (1-5) 

t~+0 7T J A 

where A is the cutoff point for which equation (jl.ip is nonoscillatory in (—00, A) and oscillatory in 
(A, 00), or equivalently, the lowest point of the essential spectrum. The spectral density function is then 
defined for A € [A, 00) by 

/(A) := p'(X) = - lira Im[m(A + ie)]. (1.6) 

7T e->0 

We now summarize some well known information on the spectrum associated with the problem 
(|1.1[) - (|1.2[1 (see, for example, [5J Chap 1-2]). For the case of periodic potentials with the above boundary 
condition, the spectrum is known to be absolutely continuous and consisting of bands interspersed with 
open intervals called gaps. For the regular periodic problem having the boundary conditions 

2/(0) = y{£), j/(0) = y'{l), (1.7) 

let the eigenvalues be ordered by Aq < Ai < • • • , where eigenvalues of multiplicity two are written twice 
in the sequence. For the regular semi-periodic problem having the boundary conditions 

2/(0) = -y(l), 2/(0) = -y'(£), (1.8) 

let the eigenvalues be ordered by /io < /ii < • • • , where eigenvalues of multiplicity two are written twice 
in the sequence. Then the eigenvalues of the periodic and semi-periodic problems occur in the order 

-00 < A = A < /io < Mi < Ai < A 2 < /j 2 < /i 3 < A 3 < A 4 • • • . 

If A denotes the self-adjoint operator associated with the SL problem (|l.l[) - (jl.2[) then the closed intervals 

[A ,/i ], [a*,Ai], [A 2 ,/j 2 ], [/i3,A 3 ],.. . (1.9) 

constitute the essential spectrum er e (or the stability set) of A. The complementary set of open intervals 

(Mo,Mi): (Ai, A 2 ), 0«2,M3), . . . (1.10) 

are the gaps in er e , or the instability set. The spectral function p(A) in (|1.5p is absolutely continuous 
and monotone increasing on o~ e , so the absolutely continuous spectrum is o~ ac (A) = a e . 

In this paper we develop new characterizations for the spectral density function /(A); one leads 
to a very efficient algorithm for its calculation. In Section 2 we summarize some general information 
concerning the SL equation (|1.1[) and a companion system of first order equations which we will utilize 
in this paper. In Section 3 the new characterizations are derived. The remaining sections develop the 
numerical scheme and show examples. 



2 Preliminaries 

In this section we give the first order system of equations which we have found to be a useful companion 
system for the study of the Sturm-Liouville equation (the PQR-cquations in (6], [7]), introduce a standard 
basis for the solution space, and state some relations which connect it to equation ([1.1)1 . 
Consider the companion first order system for U = (P, Q, R) T for A 6 (A, 00): 



dU 
dx 



d 

dx 



' p ' 




Q 




R 





A- q 
-2 2(X-q) 
0-1 





" p ' 




Q 




R 



(2.1) 



The following statements are straightforward, if occasionally tedious, to verify. 
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1. If y is any solution of the SL-equation, then ( (y 12 , —2yy, y 2 ) is a solution of equation (|2.1 

2. If we let a fundamental system of the SL-equation be defined by the initial conditions, 

" u(0,A) u(0,A) 1 _ [ 1 " 
u'(0,A) t/(0,A) J = [ 1 J ' 

then a corresponding fundamental system of solutions of equation (|2.1[) is 



(u') 2 uV {v'f 

-2uvf —[u'v + uv'] —2vv' 

„,2 „,„, „,2 



(2.2) 



(2.3) 



3. If we represent a general solution of equation (|2.1[) in the form 

a(A)E7i(K,A) + &(A)tf a (s,A) + c(A)l73(a;,A), 



P 
7? 



then (using the initial conditions (j2.2[> ) we have 

o(A) = i?(0,A), 6(A) = -Q(0, A), c(A) = P(0,A). 



(2.4) 



(2.5) 



4. The solutions {0,0} defined by the initial conditions (jl.3|) are linearly related to the solutions 
{u, v} (and vice versa) by 



9 = u cos a + v sin a, <j> — —u sin a + v cos a, 
u = 9 cos a — sin a, v — 9 sin a + </> cos a. 

5. An indefinite inner product on the solution space of equation (|2.1[) may be defined by 

(Ui, U 2 ) ■= 2(P 1 i? 2 + P2R1) — Q1Q2 = const, independent of x £ [0, 00) 

where U k = (P k , Q k , R k ), k = l,2. 

6. For any solution U — (P,Q,R) T of equation (|2.1[) , 

^(U,U) = ^[4PR-Q 2 ]=0, 

i.e. 

APR — Q 2 = const, independent of x e [0, 00) 

7. If XJ\ and U2 are any two solutions of equation (12.11) represented as in (|2.4I) . then 

(Ui,U 2 ) = 2(aic 2 + cia 2 ) - 6162 

and, in particular, 

(Ui,Ui) =4a lCl -bl 

8. If C/ = (P,Q,R) T is any solution of equation (|2.1[) expressed as in (|2.4p . then 

([/, C/) = 4Pi? - Q 2 = 4ac - b 2 . 



(2.6) 
(2.7) 



(2.8) 



(2.9) 

(2.10) 
(2.11) 

(2.12) 



9. If U = (P, Q, R) T is any solution of equation (|2.ip expressed as in (|2.4p . and if it is also written as 

U = 71 Vi + 72^2 + 73 V 3 
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where 

yf fftf W? 

-299' -\6' <j> + Oft] -2<jxf>' 

fl2 a j. j.2 



V=[V U V 2 ,V 3 ] = 



(2.13) 



is the fundamental system of (|2.1[) generated by the solutions {9, <fi} defined by (jl.3[) , then 

4ac - b 2 = 4 7 i7 3 - 7 f. (2.14) 

In fact, the result holds if 9 and are any two solutions of the SL-equation with Wa;(0(-, A), <p(-, A)) 
= 1. 

10. If y is any solution of the SL equation (|1.1[) and C/=(P,Q,R) T is any solution of companion 
system (|2.1[) then 

■fj-[Py 2 + Qyy' + R(y') 2 } =0, 

i.e., 

P(x, X)y 2 (x, A) + Q(;r, A)y(a;, X)y'(x, A) + i?(x, X)(y' (x, A) 2 = constant, independent of x. (2.15) 

In the proofs in the next sections we will make frequent use of the above results (particularly 1) 
which relate the solutions of the SL-equation to the solutions of the companion system (|2.1[) . We 
exploited similar interrelations in [6J and [7] in the study of potentials on the half line satisfying 
q G Li(0, oo). Here we make use of the same interrelations in the study of periodic potentials. 
Remark : In our previous papers [BJ, [7] the system (|2.1I) was referred to as the "PQR equations" 
(our notation); however, the analysis leading to them (particularly the motivating property (|2.15[) 
) was discovered by M. Appell (3[ in 1880. Accordingly, we will henceforth refer to this first order 
system as the Appell equations. 

3 Characterizations of the spectral density function 

In this section we give an analog of the closed form characterization obtained in [7j when q € Li(0, oo). 
For the case of a periodic potential on the half line [0, oo) the basic ideas from p], [2], [6], and [7] carry 
over, at least for values of A in the stability intervals, to yield several formulas for the spectral density 
function. 

We begin with the following definition as in pQ. 

Definition. The Sturm-Liouville equation (|1.1[) . with q(x) periodic of period £, satisfies Condition A 
for a given real value of A if and only if there exists a complex- valued solution y(x, A) for which 

y(x, A) 2 dx 

Jim = 0. (3.1) 

\y{x,\)\ 2 dx 



o 

We now have the following lemmas. 

Lemma 1. For A in the stability intervals, let 

ipi{x, A) = pi(x) exp(ik(X)x) (3-2) 

be the first Floquet solution for the characteristic exponent pi = exp(i£k(X)). Here the first Floquet 
solution for each A is understood to have the choice of k(X) such that < £k(X) < ir, and Pi(x) is 
periodic of period i. Then 

I ipi(x,X) 2 dx 

lim -\ = (3.3) 

(x, X) | dx 
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It follows from [1] Theorem 2] that the spectrum of with (| 1 . 2|) is absolutely continuous in the 

stability intervals. 

Proof. 

I ipf dx = [1 + V exp(2y£fc(A))] / exp(2ifc(A)t)p?(t) di 

o j=1 y 

1 - \exp(2i£k(X))] m f e . , / x , , 2 , , , 

1 ' exp(2ik(\)t)p 2 1 (t) dt. 



1 - exp(2itt(A)) 7 

But < fc < 7r/£, so exp(2i^fc(A)) 7^ 1 for A in the stability intervals; hence, the ratio is bounded above. 
It follows that 



r 1 



>0 



771£ — p 777 £ 

|Vi(x,A)| 2 dx / \ Pl \ 2 dt 
Jo Jo 

as m — > 00. I 

Lemma 2. Let A e U^ =0 [A2 m , /i2 m ] U [/J>2m+i,^2m+i]- Then, since Condition A holds for A in these 
stability intervals, there is a complex- valued function £(A) that is uniquely defined for A in the stability 
intervals by the properties 

(i) Im |(A) > 0, and 

(ii) 



lim 

N-yoo 



N 

{8(x,\)+£(X)(j)(x,\)) 2 dx 

~N 

\9{x,X)+£,(X)(f>(x,X)\ 2 dx 



Proof. This is proved in [TJ Lemmal]. I 

Since the solution that satisfies Condition A is unique up to a constant multiple, it follows from 
Lemma 1 and Lemma 2 that there exists a constant K 7^ such that 

ih. (x, A) = K{6(x, A) + ^X)cf>(x, A)]. (3.4) 

From (13. 2|) at x = we have 

tfji (0) cos a + ipi (0) sin a = p\ (0) cos a + [ikpi (0) + p[ (0)] sin a. 

But from (ETHl 



V>i(0) cosa + ^i(0) sina = AT[6>(0, A) + £(A)0(O, A)] cos a + K[6'(0, X) + £(A)0'(O, A)] sin a 
= X [cos 2 a — £ sin a cos a + sin 2 a + £ sin a cos a] 
= X. 

It also follows immediately from Theorem 2 of [I] that for all A in the stability intervals, we have 
that the function £(A) is, in fact, the boundary value of the Titchmarsh-Weyl to— function defined in 
tUD , that is, 

£(A) = A(X) + iB(X) := lim to(A + ie). (3.5) 
We are now ready to prove the following theorem: 

Theorem 1. For A in the stability intervals, there exists a solution U = (P, Q, R) T of Appell's equation 
(|2.ip . unique up to a constant multiple, which is periodic of period I on (0, 00). 
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Proof. For A in the stability intervals we set 



U 



' p ' 




Q 

R 





MM 



where ipi{x, A) = p\{x) exp(ifc(A)x), and ip2 = ip\- Since pi{x) is periodic of period 
is also periodic of period I. We obtain from (|3.6[) that 

K| 2 -2fc(A) ImfaS) 
-2Re(p 1 p' 1 ) 
bil 2 



" p ' 




Q 




R 





k{\?\pi 



(3.6) 



it follows that p[ 



(3.7) 



and each component is real valued and periodic with period t. 

To prove that the periodic solution is unique, up to constant multiple, consider the fundamental 
solution matrix of Appell's system of equations obtained by replacing {u,v} in (|2.3[) by the Floquet 
solutions {tpi, 1^2}, and let T(x) be the transfer matrix which carries U(x, A) to U(x + £, A), i.e. 

U(x + £, A) = T(x)U(x,X). 

Since the Floquet solutions satisfy il>j(x +£,X) = pjipj(x, A), where pj = exp(±ik(X)£), j — 1, 2, are 
the Floquet exponents for A in the stability intervals, it follows that the first and third columns of (|2.3[) 
(with tpi and ^2) are eigenvectors of T(x) with eigenvalues p\ and p\ (which are not one), and the 
second column is an eigenvector of T(x) with eigenvalue p\P2 — 1- Hence T has a one-dimensional 
eigenspace for which P, Q and R are all periodic of period £. I 

The next result provides three different representations for the spectral density function /(A). 

Theorem 2. For A in a stability interval let U — (P,Q,R) T be the periodic solution of Appell's 
system (|2.ip which is normalized by (compare (12.91) ) 

(U, U) = APR - Q 2 = 4. (3.8) 

Let {a(A), b (A) , c (A)} be the coefficients in the representation (|2.4p of this periodic solution. Then 
the spectral density function defined by (|1.6[) admits the following representations: 



/(A) = 



1 



7r[c(A) sin 2 a + 6(A) sin a cos a + a(X) cos 2 a] 

1 



7r[P(0, A) sin a — Q(0, A) sin a cos a + i?(0, A) cos a] 

1 



(3.9) 
(3.10) 
(3.11) 



Here it will be observed that the normalization (J3T8J) fixes the periodic solution only up to a ± sign; 
it is for this reason that we take the absolute value sign in these formulas to ensure that /(A) > 0, 
as required. In the applications it often happens that the denominators in the above expressions are 
positive in one stability interval and negative in another. 



tt[P(x, X)(/)(x, A) 2 + Q(x, \)<f)(x, \)(f)'(x, A) + R(x, \)<j>'(x, A) 2 



Proof. The formulas (13.91) and (|3.10l) are equivalent because the representation (|2.4j) guarantees 
that {a(A), &(A), c(A)} are given by (|2.5I) . The denominator in (|3.11l) is constant, independent of x, 
by (|2.15[) and equal to p.lOj) on evaluation at x — 0. So it suffices to prove (|3.10p subject to the 
normalization (|3.8[) . Since ipi (x, A) is linearly dependent on 9(x, A) + £(A)0(x, A) by p.4[) (where 
^(A) is the complex valued function defined on the stability intervals in Lemma 2), and ip2 = ipi is 
linearly dependent on 9(x, A) + £ (\)(j)(x, A), we may represent the periodic solution in (13.61) as 



(3.12) 





■ P(x) ' 






U(x, A) := 


Q(x) 


■= K(X) 


-[(o + + W) + (° + WW + & 




R(x) 




{0 + &){0 + &) 



G 



for some real constant K(X), independent of x. The required nomalization (|3.8[) is equivalent by 
(|2T2]) to 



4a(A)c(A) - (&(A)) 2 = 4. 



(3.13) 



Using the representation of the periodic solution in terms of the fundamental system (|2.13l) of Appell's 
equations, 

P \ 

Q =71^1+72^2+73^3, (3.14) 
Rj 

and comparing the i?-component with the i?-component in (|3.12|l . gives 

7i =K{\), 72 = 2Re(C(A))^(A), and 73 = \Z{\)\ 2 K{\). 
Hence from f|2 . 12[) and (|2.14[) we have 



APR - Q 2 = Aac - b 2 



4 7 i73 - 72 
A{K{\)) 2 [| 
4(K(A)) 2 [lm 2 £(A)] =4 



A{K{\)) 2 [|£(A)| 2 -Re 2 £(A) 



if and only if 



[X(A)Im£(A)] 2 = 1. 



(3.15) 



From (|3.5p and (jl.6p it follows (since Im£(A) > by Lemma 2) that the normalization (13.8[) holds 
if and only if 



(3.16) 



Next, we use the initial conditions (|1.3I) to evaluate the right hand side of fl3.12[) and then substitute 
into the denominator of (|3.10[) to obtain, 



7r [P(0,A) sin 2 a - Q(0, A) sin a cos a + R(0, A) cos 2 a] 

(sin 2 a + cos 2 a) ■ 1 
+ (sin 3 a cos a + sin a cos 3 a — sin 3 a cos a — sin a cos 3 a) (£ + £) 



TTi^A) 

7rif(A) 



-(2 sin 2 a cos 2 a — 2 sin 2 a cos 2 a)|£|' 
The formula (|3~TT))) now follows from f3~T^) and (|3~T7|) . ■ 



(3.17) 



A more useful characterization of /(A) is given by the following result. 
Theorem 3. Assume A is in a stability interval. Then 

VA-[u(£,X)+v'(£,X)} 2 



/(A) 



27r [u'(£, A) sin 2 a + (u(£, A) - v'(£, A)) cos a sin a - v(£, A) cos 2 a] 



(3.18) 



Here the absolute value is needed to ensure that /(A) > 0; this is due to the fact that the denominator 
could be negative in some of the stability intervals, and also corresponds to the fact that the normalization 
of {a, b, c} in (|3.19|) fixes {a, b, c} only up to a ± sign. 
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Proof. From Theorem 2 it follows that if we can construct a solution U — {P. Q, R) T of Appell's first 
order system (|2.1|) which is periodic of period I and satisfies the normalization APR — Q 2 — 4 then 
we can use it to get /(A) (e.g. from any one of the formulas (|3.9|) . (|3.10[) . or (|3.11jl ). In particular, if 
this periodic solution is represented in the form (|2.4p it follows from (|2.12[) that the coefficients {a, b, c] 
satisfy 

4a(A)c(A) - 6 2 (A) = AP{x, X)R{x, A) - Q(x, A) 2 = 4. (3.19) 

Considering only the third component of (|2.4I) it therefore suffices to generate coefficients {a, b, c} for 
which the quadratic form 

A)) 2 + bu(x, X)v(x, A) + c(v(x, A)) 2 (3.20) 

is periodic of period I and such that (|3.19[) holds for A in the stability intervals. Then /(A) is given 
by (|3.9p with this choice of {a, b, c}; or by p.lOj) . p. Ill) where {P, Q, R} is the corresponding periodic 
solution (|2.4[) of Appell's equations. To manufacture {a, b, c} consider the SL-equation (ll.lj) in the 
system form 

-(A -*(*)) l)* = A W-*> *^)=(l> 
Let ^ a (-, A) be the solution of the initial value problem 

¥' a (x) = A(x)* a (x), * a (a)=/,ae[0,oo). (3.21) 
The following facts are easily verified: 

•*>-(#)#)) <" 2 > 

[^oO^r 1 = * x (0) (3.23) 

V solutions y of (H) : ( jf^ ) = * (s) ( ^ ) (3.24) 

V solutions y o/ JED : ( ^ ) = *«(«) ( ^ ) (3.25) 

* a (as) = tf t (s) • ¥ a (t) (3.26) 
The fact that g(x) has period € implies 

*/(a; + *) = *o(aO, (3.27) 

and hence that ^ x (x + €) is periodic with period I. To generate a quadratic form in u and v which 
is periodic of period I we put 



$(x) := ( 1 )* x (x+l) ( J ) 

= (10 )*f(aj + i)* (0*x(0) 







1 jUro^foW^W" 1 ( i 
li(x) i;(x) ) 



(3.28) 



— v(x) 
u{x) 

(3.29) 

= w(l)(u(i)) 2 - \u{l) - v'(£)] u(x)v(x) - u'(£)(v(x)) 2 . (3.30) 
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Since &(x + £) — $(x) for all x £ [0, oo), we need only normalize the coefficients to achieve the 
required normalization (|3.f 9|) . Taking 7 • $>(x) so that 4ac — b 2 — 4, we find 



7 2 [-4w(*)u'(*) - - v'{t)f] (3.31) 

= 7 2 [ 4 -(«W+«'(0) 2 ] (3-32) 
= 4, (3.33) 



so that the required normalization is achieved with 



c 

and substitution of this into (|3.9|l yields (|3.18|l 



a(A) \ / -2v(*) 

6(A) = 2(u(*) -«'(*)) |, (3.34) 

(A) / V^WHAW I 2u '(£) 



4 The Numerical Method 

In this section and the following two sections we describe a new numerical algorithm for obtaining 
approximations to the spectral density function, by making use of the representation p,18|) in Theorem 3, 
and compare performance with SLEDGE. For general information and discussion of numerical methods 
for Sturm-Liouville problems we refer to Pryce's book |12j . and for the the computation of spectral 
functions using the method of SLEDGE we refer to our previous papers [TU1 HI H]- I n contrast to 
SLEDGE, the above Theorem 3 for periodic potentials enables computation of the spectral density 
function on the stability intervals by shooting (with piecewise trigonometric / hyperbolic splines) over 
a single period. 

To compute u and v we employ the method of coefficient approximation by which q is replaced by a 
step-function approximation q. We write the analog to as 

— y + q(x)y = Ay, a < x < 00 (4-1) 

and u and v will satisfy (|4.1|) with initial conditions analogous to those of u and v, respectively. We 
extend the formula (|3.18p by defining for any A 

/(A) = 2 M{0,4-^) + ^(f} (42) 

27r|u(^)sin a + (u(t) — v (£)) sin a cos a — v{£) cos a\ 

as an estimate of /. The appeal of this approach is that closed-form solutions, piecewise circular or 
hyperbolic trig functions, are known for it and v, admitting efficiencies of computation and analysis. 

There are two potential numerical challenges in trying to integrate (|4.1j) and use (|4.2j) : (1) mathe- 
matical instability when A < q(x), and (2) loss of accuracy if the numerator and demominator of (I4.2[) 
vanish simultaneously. We note that these difficulties arise in the original equations and (|3. 181) , so 
we would expect them to be inherited by any computational approach. With coefficient approximation 
it is straightforward to address both of these issues. 

In [TU] we presented a stabilizing algorithm to solve (|4.ip for the regular Sturm-Liouville problem, 
which is a boundary value problem. A similar approach will work here. First, we provide more detail of 
the algorithm. We first subdivide [0, £} into N intervals 

= X\ < x 2 < ■ ■ ■ < xn+i — t\ 

we set h n = x n+1 — x n to be the width of the nth subinterval. 

On any subinterval (x n ,x n +i) we choose q(x) = q n to be constant (usually the q value at the 
midpoint); then the differential equation (14.11) has the closed- form solution 

y{x) = y(x„)(f>' n {x - x n ) + y'(x n )(f> n (x - x n ) (4.3) 
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with 

{sinw„i/cj„ t„ > 
sinh u} n t/u n r n < 
t T n = 0, 

where 

T n = A — g„ 

and 

w„ = >/|r n |. 

It follows that 

y'(x) = -T n y(x n )(f> n (x - x n ) + y' (x n )<f>' n (x - x n ). 
In practice, one should use a truncated series expansion for small |r n |/i„, e.g., 

(f> n (t) =t[l-T n t 2 /6 + T n "t 4 /120] 

and only use the sin and sinh formulas when |r n |/i„ is sufficiently large. 
As a consequence, if we set 

Vn ■= V(x n ), Vn ■= y' ' (x n ) 



for any n, then we have the forward recurrence 



Vn+l 




Vn+l 





(j>' n {h n ) (j>n{hn) 

-T n 4> n {h n ) (j>' n {h n ) 



If we denote the coefficient matrix in (l4.6|) by A n , it is not difficult to show that it has inverse 

A~ - 



4>' n {h n ) -<f>n(hn) 
T n <Pn(h n ) 0' n (hn) 



Hence, a backward recurrence is 



Vn 




. y'n _ 





<t>n(h n ) -4>n{h n ) 
T n <t>n{h n ) (f>' n {h n ) 



Vn+1 
Vn+l 



(4.4) 



(4.5) 



(4.6) 



(4.7) 



(4.8) 



It can be seen when t„ > that A n has eigenvalues cos(u> n h n ) ± zsin(o;„/i„) and spectral radius one. 
When r n < 0, its eigenvalues are cosh(w n /i„) ± sinh(w n /i n ) and its spectral radius is exp(w n /i n ). The 
exponential factor reflects the potential mathematical instability of the initial value problem when 
A < q(x). To overcome this, define 



exp(w„/i„) r n < -e 
1 otherwise, 



and for j < k 

Introduce the scaled variables 



p(j, k) = <Jj<r j+1 



(4.9) 
(4.10) 



y n = y n /p(l,n-l) (4.11) 
y' n - y'Jp{l,n-l), (4.12) 

which satisfy the recurrences of the form (|4.6p or (I4.8[) with coefficient matrix divided by a n . These 
scaled matrices have spectral radius one. 

To use (|4.6|) requires an initial condition to start, while (|4.8I) requires a terminal condition. More 
generally, we define u F and v F to each satisfy the differential equation (14.11) with respective initial 
conditions u F (0) = 1, u F '(0) = and v F (0) = 0, v F '(0) = 1. Similarly, u B and v B satisfy the same 
differential equation but with respective terminal conditions u B (£) — 1, u B ' '(£) = and v B (£) = 0, 
v B (£) = 1. We will define the 2-vector U F to have components u F and u F ; furthermore, for n = 
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1,2,..., JV + 1 let U F denote the two-vector with components u F and u F ■ Define 2-vectors U B , V F , 
V B , U B , V F , and V B analogously. For the vectors with B superscripts, we recur backwards from 
n = N + 1 using (|4.8j) while for those with F superscripts we recur forwards from n = 1 using (|4.6[) . 
Finally, we use a tilde overscore (~) to denote the scaled versions of these recurrences. The analog of 
(|4~TT]) - (|4~T2l for Y either U or V is 

Y n F = Y F /p(l,n-l) (4.13) 
= Y B /p(n,N). (4.14) 



Since {/ and form a basis of solutions for (|4.1[) there exist constants en, C12, C21, C22 such that 

U F = Cll U B + c 12 V B (4.15) 

T/ F = c 2 iU B + C22V B (4.16) 

for every x. Define 

A = u B (v B )' - (u B )'v B , (4.17) 



then after some calculation it follows that 



Cll 


= {u F (v B y 


-(u F ) 


'v B )/A 


C12 


= (u B (u F y 


~(u B ) 


<u F )/A 


C21 


= {{v B )'v F 


- v B (v 


F y)/A 


C22 


= (u B (v F y 


~{u B ) 


'v F )/A 



(4.18) 
(4.19) 
(4.20) 
(4.21) 



for any choice of x £ [0, £}. 

From the first component of (|4.15|) 

u(£) = c llU B {£) + c 12 v B (£) = c n . 

Similarly, 



u'(£) 




f c 12 v B '(£) 


= Cl2 


v(£) 


= c 21 u B (f)4 


-C 22 V B (£) = 


: C21 




= c 21 u B '(£)- 


VC22V B \£) 


= C 2 2 


(OH 


that 







t _ ^/max{0A-[u(£)+v'{£)} 2 } 
.U A J - 



2ir\[u'(£) sin a + (u(£) - v'{£)) sin a cos a - v(£) cos a] 



\J max{0, 4 - [cu + c 22 ] 2 } 



27r|ci2sin a + (cn — C22) sin a cos a — C21 cos a| 



(4.22) 



Since the formulas (|4.18j) - (|4.21j) are valid for any x, another approach is to recur from both ends, 
computing the various u F , u B , v F , and v B , or their scaled equivalents, and 'match' at some interior 
point denoted by x = xm, to be determined below. In detail, we begin with 

uf = i, «f = o, («f)' = o, {v F y = i, 

i#+i = l> 5^+1 = 0, («W = 0, (v B +1 )' = l, 



and then compute 



^f+i = A n U F /a n (4.23) 
V;+i - A n V F /a n (4.24) 
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for n — 1, 2, . . . ,M — 1, and 



l>f = A-W^Jan (4.25) 
= A" 1 ^/^ (4.26) 
for n = N, N - 1, . . . , M. Now from (j4T3l) - (j4~T4l) and (|4~T7)) - (|4TT8)) . with n = M wc have 

rr B \l F Br F \l\ ir Br B \i r B \t B \ 
Cll = ((%) u Af - %/("m) )/(%(%) - (%) %) 

= P (i,M - i)((i&)'t& - «£(fi&)')/[p(M,iV)(u£(s£)' - 

Consequently, if we define the scale factor 

p(l,M-l) 

then en and similarly C12, C21, and C22 given in (I4.18[) - (|4.21[) must be multiplied by (m if scaled variables 
are used. To avoid rapid error buildup, it is desirable to have (m ~ 1, equivalently, 

p(l,M-l) « P (M,iV) 



« VVP) (4.28) 

with M chosen to be an index for which the approximation is best. In extreme cases, products of the 
a n may overflow, so it is best to work with their logs, i.e., from (|4.9p . h n oj n . Then (|4.28l) becomes 

M-l N 

' h n u n ^ 0.5^2' h n u n , (4.29) 

n— 1 n— 1 

where the ' on the sum means to replace h n u> n with zero for any index corresponding to r„ > 0. Care 
must also be taken in scaling en, C12, C21, and C22 by £ to do the quotient with the logs first and only 
perform the exponentiation at the end. 

Hence, for a given choice of A, the stabilized algorithm first makes an initial pass across the subin- 
tervals of [0,1] to compute the {<r„} and M. Next, the scaled forward and backward recurrences are 
performed that allow the computation of en, C12, C21, and 022- Finally (|4.22[) can be used to compute 
the estimates for /(A). This can be repeated for a sequence of ever finer meshes until convergence is 
observed. If this is not accomplished in a certain number of steps, the computation is suspended and 
an error flag is set. 

As an illustration we choose the Mathieu equation for which 

q(x) = cosx. (4.30) 

We computed the spectral function /(A) at 101 equally spaced A values in several stability intervals, with 
simple and with double shooting; the average time was measured for each method and interval. In all 
cases a hundred repetitions were made for each of the A values in order for the computer clock to produce 
a reliable time. This was done at several tolerances on the / sequence. The output is summarized in 
Table 3.1 for a Dirichlet condition at x = corresponding to the choice a — 0. Two absolute error 
tolerances were used: 10 -6 and 10 -8 . The 'failure' column shows a count of the number of values for 
which convergence was not achieved in eight mesh refinements (bisected uniform meshes). In Table 3.2 
are the corresponding values for a Neumann condition (a = 7r/2). 



Table 4.1. Simple vs. double shooting for (|4.30l) - Dirichlet. 
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simple shooting double shooting 

tolerance interval time # failures time # failures 



10~ 6 [-0.3784,-0.3476] 
[0.5949,0.9180] 
[1.2932,2.2851] 
[2.3426,4.0319] 

10~ 8 [-0.3784, -0.3476] 
[0.5949,0.9180] 
[1.2932,2.2851] 
[2.3426,4.0319] 



0.253 42 

0.080 

0.006 

0.006 

0.659 96 

0.458 51 

0.021 

0.017 



0.023 

0.007 

0.007 

0.006 

0.042 

0.013 

0.021 

0.017 



Table 4.2. Simple vs. double shooting for (|4.30j) - Neumann. 

simple shooting double shooting 

tolerance interval time # failures time # failures 



10~ 6 


[-0.3784, -0.3476] 


0.300 


47 


0.023 







[0.5949,0.9180] 


0.110 


1 


0.007 







[1.2932,2.2851] 


0.007 





0.007 







[2.3426,4.0319] 


0.009 





0.009 





io- 8 


[-0.3784, -0.3476] 


0.590 


96 


0.040 







[0.5949,0.9180] 


0.491 


68 


0.016 







[1.2932,2.2851] 


0.018 





0.018 







[2.3426,4.0319] 


0.021 





0.021 






Clearly the simple shooting approach falters on the first two stability intervals. For the other two 
intervals the quantities r n in (|4.4p are always positive so that theoretically the two methods should be 
equally reliable. The output supports this, and little time is lost from the minor overhead of the double 
shooting. For the remainder of the output in this paper the double shooting method always will be used. 

5 Indeterminate cases 

When the potential q is periodic, it is known that the spectrum exhibits spectral gaps of resolvent set 
where no spectrum can occur, i.e, where /(A) = 0. Moreover, the endpoints of a spectral gap occur at 
values of A* for which (u + v')(£, A*) = ±2, i.e., the numerator in (|3.18j) vanishes. For some examples the 
denominator may also vanish at the same A* . In such cases we might expect the numerical error to be 
large for values of A near such A*. In fact, such A* arise at endpoints of spectral gaps for any potential 
exhibiting even symmetry, i.e., q(£ — x) — q(x) for all x, as is the case for Mathieu's equation. 

Case 1 (Dirichlet): assume that for some fixed A = A* we have v(£, A*) = and u(£, A*) + v x (£, A*) = 
±2. Then, near A* we have 

4- k^,A)+^,A)] 2 = [2 + \u(£,\)+v x (t,X)\}[2-\u(t,X)+v x (£,\)\] 

= [2 + \u(£, A) + v x (£, X)\][(u x (£, A*) + v xX {£, A*))(A - A*) + 0((A - A*) 2 )] 

and 

v(£, A) = v x (£, A*)(A - A*) + 0((A - A*) 2 ), 
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so that for A < A* 



/(A) 



V[2 4 


-\u(£,\) + v x (£,X)\]\u x (£, A*)- 


f «xa(4 A*) 1 




2ir\v x {£,\*)\V\* -\ 




~ v1 2 + 


-\u(£,X) + v x (£,X)\]\u x (£, A*) - 


f «xA& A*) 1 



+ 0(X* — A) 



27r|u A (AA*)|VA^A 



(5.1) 



Case 2: Neumann: assume that for some fixed A = A* we have u x (£, A*) = and u(£, X*) + v x (£, A*) 
±2. Analogous to the Dirichlet case, we have for A > A* 



V[2 4 


- |u(4 A) + A)|]|t* A (4 A*)- 


-v xX (£, A*)| 




2ttK a (£, A*)|VA- A* 






-\u(£,X) + v x (£,X)\]\u x (£,X*)- 


rv xX (£, A*)| 


27r| Ua;A (£, A*) VA - A* 



(5.2) 

For a step-function potential the partial derivatives appearing in the above / formulas can be com- 
puted easily from the closed form solutions given in the previous section. The details are given in an 
appendix. Note that in either case we expect a 1/ •v/jA — A*] behavior near a point of indeterminacy A* . 
Our experience has shown that the expected loss of significance is not serious except (1) at very tight 
tolerances, (2) at values of A very close to A*, or (3) near gap endpoints where the gap is very narrow 
(larger A*). 

As an illustration we again choose the Mathieu potential (|4.30|) . For an absolute error tolerance of 
10~ 8 , we evaluated (|3.18l) and (|5.1I) near endpoints of the stability intervals. For (|3.18l) we also estimated 
the rate a in 

constant 
H > ~ \x- A*| Q 

by 

rater- lo sif M / f (X 1 )} 



logOAx - A*|/|A 2 - A* 



Table 5.1 displays the numerical output for a Dirichlet initial condition, where the indeterminacy occurs 
at the right-hand end of a stability interval. The respective A* values for (|5.1[) are 

{-0.347669125306, 0.918058176625, 2.28515693444}. 

Table 5.2 does the same for a Neumann initial condition, where the indeterminacy is at the left-hand 
end of a stability interval. The A* values are 

{-0.378489221265, 0.594799970122, 1.29316628334}. 
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Table 5.1. Behavior near an indeterminacy for (|4.30|) - Dirichlet. 



\ 

/\ 


f frnm \'X 1 Ah 

/ 11U111 IIO.IOIJ 




f frnm (rTTI> 

J 11 Hill \IO . Ill 


u.o^ty ( 


1 34070 
i.o^tu t y 




1 "}8fi01 


— u.o^yo 


1 ^Ofi'iO 

l.OUOOU 


u.ooi 


1 q/IRRc; 
1 .0^000 


~u.o4oy 


1 74^40 




i 78090 


-0.3485 


2.13833 


0.517 


2.16680 


-0.3481 


2.98860 


0.510 


3.00872 




1 1 OQCOf! 


0^14 


1191 8K9 
ll.ZloOZ 


Q1 ^7 


9 3^81 






Q1 fi1 


9 ^881 1 


o ^nn 


z.ooyoo 




9 001 R9 


o ^nn 


9 00981 


u.yioy 


o.oooyu 


U.oUU 


o.oo r Uo 


0.9173 


4.16048 


0.500 


4.16167 


0.9177 


6.05367 


0.500 


6.05564 


2.2831 


1.90694 




1.87531 


2.2835 


2.11789 


0.485 


2.08944 


2.2839 


2.42382 


0.488 


2.39896 


2.2843 


2.92599 


0.492 


2.90536 


2.2847 


3.99392 


0.495 


3.97862 


2.2851 


11.27750 


0.498 


11.26558 



Table 5.2. Behavior near an indeterminacy for (I4.30[) - Neumann. 



A 


/ from (13.1811 


rate 


/ from ([O]) 


-0.3784 


5.21621 


0.504 


5.22961 


-0.3780 


2.21324 


0.511 


2.23105 


-0.3476 


1.63105 


0.518 


1.65468 


-0.3472 


1.34574 


0.525 


1.37416 


-0.3468 


1.16787 


0.532 


1.20046 


-0.3464 


1.04307 




1.07943 


0.5952 


10.46971 


0.500 


10.47355 


0.5956 


7.40089 


0.501 


7.40592 


0.5960 


6.04084 


0.501 


6.04691 


0.5964 


5.22980 


0.501 


5.23678 


0.5968 


4.67613 


0.502 


4.68392 


0.5972 


4.26729 




4.27581 


1.2936 


10.78605 


0.500 


10.77975 


1.2940 


7.78143 


0.500 


7.77625 


1.2944 


6.39829 


0.499 


6.39287 


1.2948 


5.56142 


0.499 


5.55555 


1.2952 


4.98576 


0.499 


4.97941 


1.2956 


4.55873 




4.55190 
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There are slight differences between the two approaches. Since it is difficult to cafcufate exact answers 
in these cases (and A* itself), we have no easy way to judge which, if either, is more correct. By an 
inspection of intermediate quantities needed for the special formula (|5.1[) . viz., u\, v x ,x, they can be 
quite sensitive to the error in | A* — A |, as well as the tolerance. In a more positive vein, it is clear that 
the double shooting method is in agreement as to the growth rate of / near A* in the indeterminate 
situations. 



6 Other Numerical results 

In this section we exhibit computational results illustrating the algorithms developed in the previous 
sections. For brevity we choose five potentials; the first is Mathieu's equation (I4.30[) . As mentioned in 
the previous section, potentials such as this one that exhibit even symmetry have special properties. It 
can be shown that if q(i — x) = q(x) for every x, then 

u(£, A) = v'(£, A) for every A. (6.1) 

Moreover, whenever \u(£, A*)| = 1 for some A*, then either 

v(£,X*)=0 (6.2) 

or 

u'(£,X*) = 0. (6.3) 

In the case (|6.2p , A* is the left endpoint of a spectral gap when y in (|l.lj) satisfies a Dirichlet condition. 
In the case of (|6.3|) , A* is the right endpoint of a spectral gap when y in (11.11) satisfies a Neumann 
condition. 

The potentials in our other examples are 

3/(2 + sina;), (6.4) 

1/Vl - 0.75 sin 2 x, (6.5) 
(0.5 + cos x + cos2x + cos3a;)/7r, (6.6) 
sin x + 0.5 sin2x + 0.1 sin 3a:. (6-7) 

These have period £ — 2tt except for (|6.5|) that has period it. In addition to (|4.30|) . examples (|6.5|) and 
(|6.6p also have even symmetry. 

In Table 6.1a we display the endpoints of the first few stability intervals for the first two examples. 
For Mathieu's equation these are known [5] from the theory of elliptic cylinder functions. The numerical 
values agree with those found in j9], or see [H Table I]. The numerical method used was simple binary 
search (bisection) seeking the zeros of 

g(X):=2-\u(£,X)+v'(£,X)\. (6.8) 

An absolute error tolerance of 10 -8 was used in all cases. Values of / were first computed over a 
sufficiently fine grid to identify the locations of the gaps. As A increases the gap width narrows, making 
it more difficult to isolate gap boundaries. Moreover, the loss of significance in evaluating g worsens; 
eventually we may have to switch to the techniques in Section 4 to help overcome this. However, this 
was not necessary for the data in Table 5.1a. 



Table 6.1a. Stability intervals for the first two examples. 



Mathieu Example 5.4 



(-0.378489, 


-0.347669) 


( 2.250000, 


2.548882) 


( 0.594800, 


0.918058) 


( 3.055360, 


3.941647) 


( 1.293166, 


2.285157) 


( 4.146186, 


5.736211) 


( 2.342581, 


4.031922) 


( 5.796032, 


7.994726) 


( 4.035301, 


6.270837) 


( 8.010349, 


10.743819) 


( 6.270945, 


9.014297) 


(10.747778, 


13.991464) 
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Similarly, Table 6.1b contains the stability intervals for Examples (|6 . 5[) — (|6 . T|) . 

Table 6.1b. Stability intervals for the last three examples. 



Example 5.5 Example 5.6 Example 5.7 



( 1.346160, 


2.136962) 


(0.106301, 


0.247914) 


(-0.419549, 


-0.391618) 


( 2.594046, 


5.310602) 


(0.503181, 


0.995282) 


( 0.570873, 


0.840333) 


( 5.452072, 


10.356984) 


(1.311604, 


2.240365) 


( 1.362407, 


2.217768) 


(10.396276, 


17.369252) 


(2.602473, 


4.151030) 


( 2.442559, 


4.011052) 


(17.380456, 


26.372454) 


(4.198967, 


6.407883) 


( 4.078880, 


6.271355) 


(26.375745, 


37.373218) 


(6.426576, 


9.160844) 


( 6.283327, 


9.017477) 



Next we compare the new formula (|3 . 18[) with a variant of the SLEDGE code. The original SLEDGE 
[I], [8], [10] could return estimates for the spectral measure p(X) but not for the density function 
/(A); to this code was added an implementation of interpolant 3 from [11] Eqn. (3.5)], there denoted 
by (I^pi,)', in order to provide estimates for /(A). Recall that SLEDGE uses the Levitan-Levinson 
characterization of the measure p(A). This requires the calculation of many eigenvalues and suitably 
normalized eigenfunctions of a regularized Sturm-Liouville problem over a finite interval (0, b) for a 
sequence of increasingly larger b. Nevertheless, it was demonstrated in [?] that SLEDGE is capable of 
successfully handling a wide scope of problems. Our first example for Table 6.2a is the Mathieu equation 
(|4.30p with a Dirichlet initial condition (a — 0), in Table 6.2b are the data corresponding to a = 7r/6, 
and in Table 6.2c are the data corresponding to a Neumann initial condition. The internal tolerance 
used by SLEDGE in the calculation of the eigenvalues and eigenfunctions was 10 -6 , while for the new 
method it was either 10 -6 or 10 -8 , as shown. The final line of the table shows the computer time needed 
for computing / and p at 66 A points. Since SLEDGE has no choice but to compute both p and /, we 
required our new method do both as well. For brevity, only some / output values and no p values are 
shown in the tables. 

SLEDGE loses accuracy near the boundaries of the stability intervals, though it still seems to be 
converging as b increases. There is little difference in using the new formula (|3.18l) at the tighter 
tolerance, other than an increase in time. What differences there are generally occur near endpoints of 
stability intervals, especially when the formulas are indeterminate there. In all cases it is clear that the 
new approach is much faster. 
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Table 6.2a. Estimates of /(A) for Mathieu's Equation - Dirichlet. 





A 


SLEDGE 


SLEDGE 




(13.181) 




(13.181) 


-0. 


.38 





.000000 


0.000000 


0. 


.00000000 





.00000000 


-0. 


.37 





.221836 


0.221583 


0. 


.22149618 





.22149622 


-0. 


.30 





.439526 


0.438250 


0. 


.43801179 





.43801181 


-0. 


.35 


1 


.262850 


1.247271 


1 


.24515689 


1 


.24515701 


-0. 


.34 





.001181 


0.000591 





.00000000 





.00000000 










spectral g 


;ap 








0. 


.59 





.000029 


0.000006 





.00000000 





.00000000 


0. 


.GO 





.034764 


0.034936 





.03503180 





.03503178 


0. 


.70 





.176067 


0.175985 





.17595739 





.17595738 


0. 


.80 





.305406 


0.304771 





.30451657 





.30451657 


0. 


.90 





.876650 


0.855793 


o. 


.84810870 





.84810870 


0. 


.92 





.025007 


0.012603 


o. 


.00000000 





.00000000 










spectral g 


;ap 








1. 


.29 





.000896 


0.000195 


0. 


.00000000 





.00000000 


1. 


.30 





.028250 


0.036852 


0. 


.03714803 





.03714802 


1. 


.50 





.188818 


0.188736 


0. 


.18871550 





.18871550 


1. 


.75 





.268937 


0.268931 


0. 


.26892936 





.26892936 


2. 


.00 





.336884 


0.336798 





.33675478 





.33675478 


2. 


.25 





.673109 


0.591205 


0. 


.56713172 





.56713172 


2. 


.28 





.391315 


0.839088 


1 


.23367035 


1 


.23367034 


2. 


.30 





.267555 


0.179038 


o. 


.00000000 





.00000000 










spectral g 


;ap 








2. 


.34 





.184005 


0.071690 


0. 


.00000000 





.00000000 


2. 


.50 





.330775 


0.330639 


0. 


.33062792 





.33062792 


2. 


.75 





.392572 


0.392579 


0. 


.39258965 





.39258966 


3. 


.00 





.431331 


0.431330 





.43133033 





.43133034 


3. 


.25 





.463740 


0.463753 





.46375540 





.46375540 


3 


.50 





.493118 


0.493118 





.49311812 





.49311813 




b 




64tt 


1287T 










tolerance 




10~ 6 


10~ 6 




10~ 6 




io- 8 


total time 




72.67 


246.35 




0.10 




0.17 
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Table 6.2b. Estimates of /(A) for Mathieu's Equation - a = tt/6. 





A 


SLEDGE 


SLEDGE 




(|3.18D 




(13.181) 


-0. 


.38 





.000000 


0.000000 


0. 


.00000000 





.00000000 


-0. 


.37 





.254422 


0.254359 





.25428685 





.25428585 


-0. 


.36 





.358362 


0.358111 





.35803367 





.35803367 


-0. 


.35 





.271701 


0.271974 





.27213584 





.27213584 


-0. 


.34 





.000000 


0.177284 


0. 


.00000000 





.00000000 










spectral g 


ap 








0. 


.59 





.000000 


0.001934 





.00000000 





.00000000 


0. 


.GO 





.046178 


0.046391 





.04652124 





.04652122 


0. 


.70 





.213029 


0.212952 


0. 


.21292212 





.21292210 


0. 


.80 





.311584 


0.311230 





.31111121 





.31111121 


0. 


.90 





.334158 


0.336703 





.33591487 





.33591485 


0. 


.92 





.064154 


0.065358 





.00000000 





.00000000 










spectral g 


ap 








1. 


.29 





.018331 


0.015724 


0. 


.00000000 





.00000000 


1. 


.30 





.036250 


0.048914 





.04930685 





.04930685 


1. 


.50 





.225283 


0.225239 


0. 


.22523166 





.22523166 


1. 


.75 





.289651 


0.289648 


0. 


.28965416 





.28965416 


2. 


.00 





.327035 


0.327015 


0. 


.32700588 





.32700588 


2. 


.25 





.365013 


0.371550 


0. 


.36740588 





.36740588 


2. 


.28 





.550911 


0.173426 


0. 


.27382995 





.27382995 


2. 


.30 


1 


.003042 


0.719051 





.00000000 





.00000000 










spectral g 


ap 








2. 


.34 





.604674 


0.435577 


0. 


.00000000 





.00000000 


2. 


.50 





.324044 


0.324191 





.32423291 





.32423291 


2. 


.75 





.347326 


0.347336 





.34733465 





.34733465 


3 


.00 





.356739 


0.356748 


0. 


.35675152 





.35675152 


3 


.25 





.362111 


0.362122 





.36212172 





.36212172 


3 


.50 





.365271 


0.365276 





.36527626 





.36527626 




b 




64tt 


1287T 










tolerance 




10~ 6 


10~ 6 




10~ 6 




io- 8 


total time 




81.45 


231.55 




0.09 




0.14 
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Table 6.2c. Estimates of /(A) for Mathieu's Equation - Neumann. 





A 


SLEDGE 


SLEDGE 




(|3.18D 




(13.181) 


-0 


38 





000000 


0.000000 





00000000 





00000000 


-0 


37 





458625 


0.457836 





45743969 





45743978 


-0 


36 





231738 


0.231458 





23132070 





23132067 


-0 


35 





081480 


0.081348 





08137240 





08137222 


-0 


34 





000162 


0.000080 





00000000 





00000000 










spectral g 


ap 











59 





019597 


0.009810 





00000000 





00000000 





60 


2 


988328 


2.893315 


2 


89226446 


2 


89226447 





70 





577118 


0.576146 





57582799 





57582799 





80 





333326 


0.332866 





33272798 





33272798 





90 





119832 


0.119566 





11946722 





11946721 





92 





000890 


0.003814 





00000000 





00000000 










spectral g 


ap 








1 


29 





113623 


0.057586 





00000000 





00000000 


1 


30 


2 


917653 


2.769684 


2 


72749873 


2 


72749865 


1 


50 





538198 


0.537216 





53689910 





53689910 


1 


75 





376904 


0.376796 





37675762 





37675762 


2 


00 





300894 


0.300884 





30087526 





30087526 


2 


25 





162234 


0.180877 





17865547 





17865547 


2 


28 





062850 


0.040850 





08212985 





08212987 


2 


30 





160618 


0.097636 





00000000 





00000000 










spectral g 


ap 








2 


34 





356154 


0.299811 





00000000 





00000000 


2 


35 





405038 


0.440822 





82126926 





82126926 


2 


50 





307667 


0.306737 





30645078 





30645078 


2 


75 





258179 


0.258104 





25808419 





25808419 


3 


00 





234940 


0.234913 





23490391 





23490391 


3 


25 





218489 


0.218484 





21847979 





21847979 


3 


50 





205484 


0.205474 





20547041 





20547041 




b 




64tt 


1287T 










tolerance 




10~ 6 


10~ 6 




10~ 6 




io- 8 


total time 




73.64 


239.55 




0.09 




0.16 



7 Appendix: estimating variational quantities 



Here we derive a method for computing the partial derivative with respect to A of the quantities given 
in Section 5 for overcoming indeterminacies. Recall from Section 4 that the forward recurrence is 



u n+l — StnUn 



so that 



Omitting the n subscripts for now, we have 



(7.1) 
(7.2) 



A, 



<Px\ <P\ 

-t x 4> - t4>\ 4> x \ 



But from (gl} 

r A = 1; 

furthermore, for either sign on r it is easily shown that at x = x n+ i 

4>x\ = -h(f>/2, 
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and 

4>\ = {H x -4>)/{2t). 

Consequently, with the n subscripts restored, a forward recursion is 

(h n ) (h n <j> n , x (h n ) - <j> n (h n ))/T n 

-<f>n(h n ) - h„cf> ntX (hn) -h n (j) n (h n ) 

(h n ) <t>n{h n ) 1 dU n 

-T n (j>n(hn) 4>n.x(h n ) 8X 



dX 



= 0.5 



U, 



(7.3) 



with 



dUf 
dX 



The forward recurrence for V F X is identical - only the initial conditions on V F differ. 
A similar analysis using the A~ l as given in (|4.7[) leads to the backward recurrence 



dul 

dX 



= 0.5 



(h n ) (4>n(h n ) - hn<j> n ,x(h n ))/r n 

4*n(hn) + h n (f>n }X (h n ) (K) 

4>n,x(h n ) ~(t>n(hn) 
Tn4>n(h n ) 4>n,x(h n ) 



u 



n+1 



dX 



(7.4) 



with 



dX 



Again this holds as well for V B X with appropriate terminal values for V B +1 . 

As we saw in Section 4 it is desirable to scale the variables. We will use the same notation as §4, 
and will first develop the formulas for U F as those for U B , V F , and V B are analogous. Following 
we define 

= U F /p(l,n-l) 
Un = U, B /p(n : N). 

Differentiation with respect to A yields 



But 



dX dX 
dp(l,n-l) 



" /p(l, n-l)-U£ '-/p(l, n - iy. 



dX 



ox 



/p(l,n-iy 



[o.5Mi,«-i) 2 lE 



hjp(l,n - 1) 



= ]>>,/v/^)/(2p(l,n-l)), 
where the sum is taken over j, 1 < j < n, for which Tj < 0. Hence, 



and similarly 



dx 

dV F 
dX 

dX 

dv B 

dX 



1 



p(l, 71—1) 
1 

P(l,n-1) 
1 



9UF_ F \ - h 3 



9A 



0V F 

^ v n 

dX 



p(n, AT) 
1 

P(",A0 



<9A 
<9A 



(7.5) 

(7.6) 

(7.7) 

(7.8) 
(7.9) 
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The sums for the -F-superscripted cases run from j = ltoj = n—l, while those for the ^-superscripted 
cases go from j = n to j — N. In either case, the indices for which tj > are omitted. Recall that we 
expect a more stable algorithm if we recur with the scaled variables, for example, using 



dX 



= (0.5/<r„) 

+(iK) 



(h n ) (h n <f> n ,x(h n ) - 4>n{h n ))/T n 

-<t>n.{h n ) - h n <t>n,x{hn) -h n (j) n {h n ) 



U, 



{h n ) <t>n{h n ) 
-T n (t)n{h n ) 4> n ,x(h n ) 



dul 

dX 



instead of (|7.3p . For the forward recurrences we take n = 1,2, 
recurrences n = N, N — 1, . . . , M. 

It remains to recover the desired values 



u x (£) 


dU£ +1 


u x ,\{£) 


dX 


vx(£) 




v x ,x(£) _ 


dX 



from the scaled variables. 

From the first component of f|4. 1 5[) we have for any x £ [0, £} 



Ux 



Similarly 



b i b , B® c n 
cnu A +ci 2 v x +u 



den 
dX 



at x = 



Ux\ = 



B B B 

CnU xX + c 12 v xX + u x 



,M 



(7.10) 

1, while for the backward 



dX 



dc 



12 



dX 



at x = 



; dc 12 
dX 



,dci2 
dX 



and 



Finally 



and 



Similarly, 



B | B . B^ c 21 . B<9c 2 2 

V\ = C 2 1U X + C 22 V X + U + V 

dc 21 



dX 



<x\> x — £. 



v x x = c 21 u xX + c 22 v xX + u x — + V x — 
dc 22 



dX 



at x = 



ox y x x ' dx 

OA 



_ B„,B i „.B„,B „,B „,B „,B„,B 

dX 



= uZvZ + u"vZ x -uZ x v°-uZvi 



^ = - «f- F )^/A 2 + [uM + u B u F xX - u B xX u F - „M ]/A 
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and 

<9C22 i B F B F\ ® / A 2 . r B F . B F B F B Fi / a 

"aF = _ ( u ^ )gX/ A « x a- w xA« -u x v x \/A. 

These last five are all to be evaluated at x = Xm ■ By inspecting the scale factors, it follows that, as in 
Section 2, we must multiply by (m given by (|4.27p when using the scaled variables. 

After the double-shooting and matching with scaled variables, the results are to be substituted into 
(|5.ip or (|5.2p . While these formulas seem complicated, for this paper they are only to be used in the 
neighborhood of a 0/0. The computation of A* itself may be done using a characterization in terms of 
eigenvalues (see [3]), or by searching for zeros of the numerators in the expressions for /, (|3.18|) . i.e., 
zeros of 2 - \u F {l, A) + v F '(£, A)| = 2 - |c u + c 22 |. 

Since we have no test problems with closed form solutions to verify computer output for the vari- 
ational variables, we have compared our algorithm with finite difference approximations. Table 7.1 
contains data for several of our examples with a Dirichlet initial condition. In all cases a central differ- 
ence was used with a stepsize of 10 , and an absolute error tolerance of 10~ 8 was used for u x \ and ua- 
The agreement is good in all cases. 



Table 7.1. Finite difference estimates compared to u x \ and v\. 



Example 


A 


A\U X 


u x \ 




v\ 


(|4.30p 


-0.35 


-63.7915 


-63.7916 


-56.2402 


-56.24019 




1.00 


-1.6844 


-1.684311 


-5.2775 


5.277455 




2.00 


1.3092 


1.309169 


-2.2702 


-2.270148 




2.00 


-1.7131 


-1.713098 


2.3125 


2.312439 




3.00 


0.9515 


0.9514705 


0.0938 


0.009380 




5.00 


-2.6099 


-2.609927 


0.6906 


0.690553 




-0.40 


-5.870 


-5.870013 


-112.35 


-112.3457 




1.00 


2.6079 


-2.607899 


5.5210 


5.521029 




2.00 


1.9781 


1.978065 


-1.8360 


-1.836113 



We conclude this section with numerical data illustrating the overhead required for the additional 
calculation of the variational variables u F , u F x , v F and v F x . Table 7.2 has timing data for all five 
examples where 601 /(A) evaluations were made over a uniform grid of A in the intervals shown. Examples 
(|4.30[) . (|6.5[) . and (|6.7|) had Dirichlet initial conditions; the other two had Neumann. The column labelled 
'basic' gives the time required for just the u F , u x , v F , and v x calculations; the final column gives the 
time required to compute all eight variables. 
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Table 7.2. Timings for calculation of basic and variational solutions. 



Example 


Interval 


Tolerance 


Basic 


All 


(|4.aop 


[1.2, 7.2] 


10" 4 


0.027 


0.031 






10" 6 


0.048 


0.086 






10" 8 


0.103 


0.120 






io- 10 


0.126 


0.149 


EH 


[3.0, 9.0] 


10" 4 


0.036 


0.041 






10~ 6 


0.078 


0.121 






io- 8 


0.146 


0.164 






ht 10 


0.180 


0.203 




[2.0, 8.0] 


io- 4 


0.040 


0.046 






10~ 6 


0.090 


0.117 






io- 8 


0.162 


0.178 






HT 10 


0.198 


0.222 




[1.2, 7.2] 


io- 4 


0.048 


0.054 






10~ 6 


0.143 


0.162 






io- 8 


0.196 


0.214 






io- 10 


0.368 


0.524 




[1.5, 7.5] 


IO" 4 


0.050 


0.056 






10~ 6 


0.138 


0.169 






io- 8 


0.204 


0.224 






io- 10 


0.336 


0.510 


Totals 




IO" 4 


0.201 


0.228 






IO" 6 


0.497 


0.655 






10~ 8 


0.811 


0.900 






io- 10 


1.208 


1.608 



Despite the doubling in the number of dependent variables, the overhead is increased by only 10% 
to 30% in the totals, largely because the basic and variational variables use the same transcendental 
function values for a fixed A. 
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